Critical behavior and Griffiths effects in the disordered contact process 
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We study the nonequilibrium phase transition in the one-dimensional contact process with 
quenched spatial disorder by means of large-scale Monte-Carlo simulations for times up to 10 9 and 
system sizes up to 10 7 sites. In agreement with recent predictions of an infinite-randomness fixed 
point, our simulations demonstrate activated (exponential) dynamical scaling at the critical point. 
The critical behavior turns out to be universal, even for weak disorder. However, the approach to 
this asymptotic behavior is extremely slow, with crossover times of the order of 10 4 or larger. In 
the Griffiths region between the clean and the dirty critical points, we find power-law dynamical 
behavior with continuously varying exponents. We discuss the generality of our findings and relate 
them to a broader theory of rare region effects at phase transitions with quenched disorder. 
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I. INTRODUCTION 

The nonequilibrium behavior of many-particle systems 
has attracted considerable attention in recent years. Of 
particular interest are continuous phase transitions be- 
tween different nonequilibrium states. These transitions 
are characterized by large scale fluctuations and collec- 
tive behavior over large distances and times very similar 
to the behavior at equilibrium critical points. Exam- 
ples of such nonequilibrium transitions can be found in 
population dynamics and epidemics, chemical reactions, 
growing surfaces, and in granular flow and traffic jams 
(for recent reviews see, e.g., Refs. HHHSHBQ) 

A prominent class of nonequilibrium phase transitions 
separates active fluctuating states from inactive, absorb- 
ing states where fluctuations cease entirely. Recently, 
much effort has been devoted to classifying possible uni- 
versality classes of these absorbing state phase transitions 
0, @ . The generic universality class is directed perco- 
lation (DP) y|. According to a conjecture by Janssen 
and Grassberger 9], all absorbing state transitions with 
a scalar order parameter, short-range interactions, and 
no extra symmetries or conservation laws belong to this 
class. Examples include the transitions in the contact 
process 0, catalytic reactions ^lj, interface growth 
|12| , or turbulence [r| . In the presence of conservation 
laws or additional symmetries, other universality classes 
can occur, e.g., the parity conserving class [ill ITU Il6| 
or the ^-symmetric directed percolation (DP2) class 

In realistic systems, one can expect impurities and de- 
fects, i.e., quenched spatial disorder, to play an impor- 
tant role. Indeed, it has been suggested that disorder 
may be one of the reasons for the surprising rarity of ex- 
perimental realizations of the ubiquitous directed perco- 
lation universality class |2(| . The only verification so far 
seems to be found in the spatio-temporal intermittency 
in ferrofluidic spikes [2l|. 

The investigation of disorder effects on the DP transi- 
tion has a long history, but a coherent picture has been 
slow to emerge. According to the Harris criterion |22T . l23| . 



a clean critical point is stable against weak disorder if the 
spatial correlation length critical exponent v±_ fulfills the 
inequality 



dv±_ > 2, 



(1) 



where d is the spatial dimensionality. The DP univer- 
sality class violates the Harris criterion in all dimensions 
d < 4, because the exponent values are v± ~ 1.097 (ID), 
0.73 (2D), and 0.58 (3D) 0. A field-theoretic renor- 
malization group study [24| confirmed the instability of 
the DP critical fixed point. Moreover, no new critical 
fixed point was found. Instead the renormalization group 
displays runaway flow towards large disorder, indicat- 
ing unconventional behavior. Early Monte-Carlo simu- 
lations 23] showed significant changes in the critical ex- 
ponents while later studies [25| of the two-dimensional 
contact process with dilution found logarithmically slow 
dynamics in violation of power-law scaling. In addition, 
rare region effects similar to Griffiths singularities j2|J 
were found to lead to slow dynamics in a whole pa- 
rameter re gion in the vicinity of the phase transition 

mm mm up. 

Recently, an important step towards understanding 
spatial disorder effects on the DP transition has been 
made by Hooyberghs et al. [3l| . These authors used the 
Hamiltonian formalism [3^ to map the one-dimensional 
disordered contact process onto a random quantum 
spin chain. Applying a version of the Ma-Dasgupta- 
Hu strong-disorder renormalization group 33], they 
showed that the transition is controlled by an infinite- 
randomness critical point, at least for sufficiently strong 
disorder. This type of fixed point leads to activated 
(exponential) rather than power-law dynamical scaling. 
For weaker disorder, Hooyberghs et al. 0| used com- 
puter simulations and predicted non-universal continu- 
ously varying exponents, with either power-law or expo- 
nential dynamical correlations. 

In this paper, we present the results of large-scale 
Monte-Carlo simulations of the one-dimensional contact 
process with quenched spatial disorder. Using large sys- 
tems of up to I0 7 sites and very long times (up to I0 9 ) 
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we show that the critical behavior at the nonequilib- 
rium phase transition is indeed described by an infinite- 
randomness fixed point with activated scaling. More- 
over, we provide evidence that this behavior is universal, 
i.e., it occurs even in the weak-disorder case. However, 
the approach to this universal asymptotic behavior is ex- 
tremely slow, with crossover times of the order of 10 or 
larger which may explain why nonuniversal (effective) ex- 
ponents have been seen in previous work. We also study 
the Griffiths region between the clean and the dirty crit- 
ical points. Here, we find power-law dynamical behavior 
with continuously varying exponents, in agreement with 
theoretical predictions. 

This paper is organized as follows. In section [HJ we in- 
troduce the model. We then contrast power-law scaling 
as found at conventional critical points with activated 
scaling arising from infinite-randomness critical points. 
We also summarize the predictions for the Griffiths re- 
gion. In section 11111 we present our simulation method 
and the numerical results together with a comparison 
to theory. We conclude in section IIVI by discussing the 
generality of our findings and their relation to a broader 
theory of rare region effects at phase transitions with 
quenched disorder. 



II. THEORY 
A. Contact process with quenched spatial disorder 

We start from the clean contact process ^(j; a P ro ~ 
totypical system in the DP universality class. It can be 
interpreted, e.g., as a model for the spreading of a dis- 
ease. The contact process is defined on a d-dimensional 
hypercubic lattice. Each lattice site r can be active (oc- 
cupied by a particle) or inactive (empty) . In the course of 
the time evolution, active sites can infect their neighbors, 
or they can spontaneously become inactive. Specifically, 
the dynamics is given by a continuous-time Markov pro- 
cess during which particles are created at empty sites at 
a rate An/ (2c?) where n is the number of active nearest 
neighbor sites. Particles are annihilated at rate /j, (which 
is often set to unity without loss of generality). The ratio 
of the two rates controls the behavior of the system. 

For small birth rate A, annihilation dominates, and the 
absorbing state without any particles is the only steady 
state (inactive phase). For large birth rate A, there is a 
steady state with finite particle density (active phase). 
The two phases are separated by a nonequilibrium phase 
transition in the DP universality class at A = A°. The 
central quantity in the contact process is the average den- 
sity of active sites at time t 

<>(*) = £3 £M*)> ( 2 ) 

r 

where n r (t) is the particle number at site r and time t, 
L is the linear system size, and (. . .) denotes the average 



over all realizations of the Markov process. The longtime 
limit of this density (i.e., the steady state density) 

p s tat = hm p(t) (3) 

t — >oo 

is the order parameter of the nonequilibrium phase tran- 
sition. 

Quenched spatial disorder can be introduced by mak- 
ing the birth rate A a random function of the lattice site 
r. We assume the disorder to be spatially uncorrelated; 
and we use a binary probability distribution 

P[A(r)] = (1 - p) 8[X(r) - A] + pS[X(r) - cA] (4) 

where p and c are constants between and 1. This dis- 
tribution allows us to independently vary spatial density 
p of the impurities and their relative strength c. The 
impurities locally reduce the birth rate, therefore, the 
nonequilibrium transition will occur at a value A c that is 
larger than the clean critical birth rate A°. 

B. Conventional power-law scaling 

In this subsection we summarize the phenomenological 
scaling theory of an absorbing state phase transition that 
is controlled by a conventional fixed point with power-law 
scaling (see, e.g., Ref. jf|). The clean contact process falls 
into this class. 

In the active phase and close to critical point A c , the 
order parameter p sta t varies according to the power law 

Art* ~ (A — X c f ~ A* 3 (5) 

where A = (A — A c )/A c is the dimensionless distance 
from the critical point, and is the critical exponent of 
the particle density. In addition to the average density, 
we also need to characterize the length and time scales 
of the density fluctuations. Close to the transition, the 
correlation length £j_ diverges as 

a ~ |A|-"- . (6) 

The correlation time £m behaves like a power of the cor- 
relation length, 

£||~£L (?) 

i.e., the dynamical scaling is of power-law form. Conse- 
quently the scaling form of the density as a function of 
A, the time t and the linear system size L reads 

p(A, t, L) = 6 /3/ ^p(A6~ 1/ ^ , tb\ Lb) . (8) 

Here, b is an arbitrary dimensionless scaling factor. 

Two important quantities arise from initial conditions 
consisting of a single active site in an otherwise empty 
lattice. The survival probability P s describes the proba- 
bility that an active cluster survives when starting from 
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such a single-site seed. For directed percolation, the sur- 
vival probability scales exactly like the density [34| . 

P s (A,t,L) = b^/ v ^P s {Ab- 1 '^,tb z ,Lb) . (9) 

Thus, for directed percolation, the three critical expo- 
nents /3, 2/j_ and z completely characterize the critical 
point. 

The pair connectedness function C(r',t',r,t) = 
{n T '(t') n r (t)) describes the probability that site r' is ac- 
tive at time t' when starting from an initial condition 
with a single active site at r and time t. For a clean sys- 
tem, the pair connectedness is translationally invariant 
in space and time. Thus, it only depends on two argu- 
ments C(r, t', r, t) = C(r' — r, t — t'). Because C involves 
a product of two densities, its scale dimension is 2(3 /i/±, 
and the full scaling form reads [3j| 

C(A, r, t, L) = b 2p/v ^C{Ab- 1/u ^- , rb, tb z , Lb) . (10) 

The total number of particles N when starting from a 
single seed site can be obtained by integrating the pair 
connectedness C over all space. This leads to the scaling 
form 

N(A, t, L) = b 2f3 / u± - d N(Ab- 1/v± , tb z , Lb) . (11) 

At the critical point, A — 0, and in the thermody- 
namic limit, L — > oo, the above scaling relations lead 
to the following predictions for the time dependencies 
of observables: The density and the survival probability 
asymptotically decay like 

p{t)~t~ s , P s (t)~r 5 (12) 

with S = f3/{y±z). In contrast, the number of particles 
in a cluster starting from a single seed site increases like 

N(t) ~ t e (13) 

where = d/z — 2f3/(v±z) is the so-called critical initial 
slip exponent. 

Highly precise estimates of the critical exponents for 
clean one-dimensional directed percolation have been ob- 
tained by series expansions |3q|: /3 = 0.276486, v±_ = 
1.096854, z = 1.580745, S = 0.159464, and 6 = 0.313686. 

C. Activated scaling 

In this subsection we summarize the scaling theory for 
an infinite-randomness fixed point with activated scaling, 
as has been predicted to occur in absorbing state tran- 
sitions with quenched disorder j3l]]. It is similar to the 
scaling theory for the quantum phase transition in the 
random transverse field Ising model |37| . 

At an infinite-randomness fixed point, the dynamics is 
extremely slow. The power-law scaling Q gets replaced 
by activated dynamical scaling 

M£||)~& (14) 



characterized by a new exponent ip. This exponential 
relation between time and length scales implies that the 
dynamical exponent z is formally infinite. In contrast, 
the static scaling behavior remains of power law type. 

Moreover, at an infinite-randomness fixed point the 
probability distributions of observables become ex- 
tremely broad, so that averages are dominated by rare 
events such as rare spatial regions with large infection 
rate. In such a situation, averages and typical values of a 
quantity do not necessarily agree. Nonetheless, the scal- 
ing form of the average density at an infinite-randomness 
critical point is obtained by simply replacing the power- 
law scaling combination tb z by the activated combination 
\n(t)b^' in the argument of the scaling function: 

p(A,ln(t),i) = b /3 / Ux p(Ab- 1 / Ux ,lii{t)b^,Lb) . (15) 

Analogously, the scaling forms of the average survival 
probability and the average number of sites in a cluster 
starting from a single site are 

P s (A,ln(t),L) = 6^P s (A6- 1 /^,l n (t)6^,L&) (16) 
N(A,ln(t),L) = b 20 ' v ^- d N{Ab- x/v ^ M{t)b^ ,Lb\YI) 

These activated scaling forms lead to logarithmic time 
dependencies at the critical point (in the thermodynamic 
limit). The average density and the survival probability 
asymptotically decay like 

p(t) ~ lHt)]- S , P s (t) ~ [\n(t)}- S (18) 

with S = f3/(i/±ip) while the average number of particles 
in a cluster starting from a single seed site increases like 

N(t) ~ lln(t)f (19) 

with = d/ijj — 2P/(u±ip). These are the relations we 
are going to test in this paper. 

Within the strong disorder renormalization group ap- 
proach of Ref. [3l|, the critical exponents of the disor- 
dered one-dimensional contact process can be calculated 
exactly. Their numerical values are (3 = 0.38197, v±_ = 2, 
^ = 0.5, 5 = 0.38197, and = 1.2360. 

D. Griffiths region 

The inactive phase of our disordered contact process 
with the impurity distribution Q can be divided into 
two regions. For birth rates below the clean critical 
point, A < A°, the behavior is conventional. The sys- 
tem approaches the absorbing state exponentially fast in 
time. The decay time increases with A and diverges as 
| A — A°|~ ziyi where z and i/± are the exponents of the 
clean critical point [H, |3!| . 

The more interesting region is the so-called Griffiths 
region j2^, H(j which occurs for birth rates between the 
clean and the dirty critical points, A" < A < A c . The sys- 
tem is globally still in the inactive phase, i.e., the system 
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eventually decays into the absorbing state. However, in 
the thermodynamic limit, one can find arbitrarily large 
spatial regions devoid of impurities. For A° < A < A c , 
these so-called rare regions are locally in the active phase. 
Because they are of finite size, they cannot support a non- 
zero steady state density but their decay is very slow be- 
cause it requires a rare, exceptionally large density fluc- 
tuation. 

The contribution of the rare regions to the time evolu- 
tion of the density can be estimated as follows [2^, [27J • 
The probability w for finding a rare region of linear size 
L r devoid of impurities is (up to pre-exponcntial factors) 
given by 

w{L r ) ~ exp(-pL*) (20) 

where p is a nonuniversal constant which for our binary 
disorder distribution is given by p = — ln(l — p). The 
long-time decay of the density is dominated by these rare 
regions. To exponential accuracy, the rare region contri- 
bution to the density can be written as 

p(t) ~ J dL r L d r w(L r )exp [-t/r(L r )] (21) 

where r(L r ) is the decay time of a rare region of size 
L r . Let us first discuss the behavior at the clean critical 
point, A 9 , i.e., at the boundary between the conventional 
inactive phase and the Griffiths region. At this point, 
the decay time of a single, impurity-free rare region of 
size L r scales as r(L r ) ~ L z r as follows from finite size 
scaling |4jj. Here z is the clean critical exponent. Using 
the saddle point method to evaluate the integral 1)21 [I. 
we find the leading long-time decay of the density to be 
given by a stretched exponential, 

lnp(t) ~ -p z/(d+z) t d/{d+z) , (22) 

rather than a simple exponential decay as for A < \ a c . 

Inside the Griffiths region, i.e., for A° < A < A c , the 
decay time of a single rare region depends exponentially 
on its volume, 

T(L r ) ~ exp(aLf) (23) 

because a coordinated fluctuation of the entire rare re gion 
is required to take it to the absorbing state |23l 1271 |41| . 
The nonuniversal prefactor a vanishes at the clean critical 
point A° and increases with A. Close to A°, it behaves 
as a ~ £j d ~ (A — A") dl/J - with v\_ the clean critical 
exponent. Repeating the saddle point analysis of the 
integral (|21|l for this case, we obtain a power-law decay 
of the density 

p(t) ~ t-p' a = r d ' z ' (24) 

where z' = da/p is a customarily used nonuniversal dy- 
namical exponent in the Griffiths region. Its behav- 
ior close to the dirty critical point A c can be obtained 
within the strong disorder rcnormalization group method 
[3TI . When approaching the phase transition, z' di- 
verges as z' ~ I A — \ c \~^ u± - where tp and v± are the 
exponents of the dirty critical point. 




P c 



FIG. 1: Phase diagrams of the disordered contact process. 
Left: Birth rate A vs. impurity concentration p for fixed im- 
purity strength c = 0.2. Right: A vs. c for fixed p — 0.3. 



III. MONTE-CARLO SIMULATIONS 
A. Method and overview 

We now turn to the main part of the paper, extensive 
Monte-Carlo simulations of the one-dimensional contact 
process with quenched spatial disorder. There is a num- 
ber of different ways to actually implement the contact 
process on the computer (all equivalent with respect to 
the universal behavior). We follow the widely used al- 
gorithm described, e.g., by Dickman jl^. Runs start 
at time t = from some configuration of occupied and 
empty sites. Each event consists of randomly selecting an 
occupied site r from a list of all N p occupied sites, select- 
ing a process: creation with probability A(r)/[1 + A(r)] 
or annihilation with probability 1/[1 + A(r)] and, for cre- 
ation, selecting one of the neighboring sites of r. The 
creation succeeds, if this neighbor is empty. The time 
increment associated with this event is 1/N p . Note that 
in this implementation of the disordered contact process 
both the creation rate and the annihilation rate vary from 
site to site in such a way that their sum is constant (and 
equal to one). 

Using this algorithm, we have performed simulations 
for system sizes between L — 1000 and L = 10 7 . We have 
studied impurity concentrations p = 0.2,0.3,0.4,0.5,0.6 
and 0.7 as well as relative impurity strengths of c = 
0.2,0.4,0.6 and 0.8. To explore the extremely slow dy- 
namics associated with the predicted infinite-randomness 
critical point, we have simulated very long times up to 
t = 10 9 which is, to the best of our knowledge, at least 
three orders of magnitude in t longer than previous sim- 
ulations of the disordered contact process. In all cases 
we have averaged over a large number of different disor- 
der realizations, details will be mentioned below for each 
specific set of calculations. 

Figure Ogives an overview over the phase diagram re- 
sulting from our simulations. As expected, the critical 
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FIG. 2: Overview of the time evolution of the density for a 
system of 10 6 sites with p = 0.3 and c = 0.2. The clean critical 
point A° w 3.298 and the dirty critical point A c ~ 5.24 are 
specially marked. 
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FIG. 3: Time evolution of the density at the clean critical 
point = 3.298 for systems of 10 7 sites with c = 0.2 and 
several p. The straight lines are fits to the stretched exponen- 
tial lnp(t) ~ -Et d/{d+z) predicted in eq. witn d = 1 and 
the clean z = 1.580. Inset: Decay constant E vs. p z ^ d+z \ 



birthrate A c increases with increasing impurity concen- 
tration p. It also increases with decreasing birth rate 
on the impurities, i.e., a decreasing relative strength c. 
For p = or c = 1, we reproduce the well-known clean 
critical birth rate X° c ps 3.298 HJ. In the following sub- 
sections we discuss the behavior in the vicinity of the 
phase transition in more detail. 



B. Time evolution starting from full lattice 

In this subsection we discuss simulations which follow 
the time evolution of the average density starting from a 
full lattice. This means, at time t = 0, all sites are active 
and p(0) = 1. 

Figure |2 gives an overview of the time evolution of the 
density for a system of 10 6 sites with p = 0.3,c=0.2, cov- 
ering the A range from the conventional inactive phase, 
A < A^ all the way to the active phase, A > A c . The data 
are averages over 480 runs, each with a different disorder 
realization. For birth rates below and at the clean critical 
point A° ~ 3.298, the density decay is very fast, clearly 
faster then a power law. Above A°, the decay becomes 
slower and asymptotically seems to follow a power-law. 
For even larger birth rates the decay seems to be slower 
than a power law while the largest birth rates give rise to 
a nonzero steady state density, i.e., the system is in the 
active phase. 



1. Griffiths region 

Let us investigate the different parameter regions in 
more detail, beginning with the behavior at the clean crit- 
ical point A°, i.e., at the boundary between the Griffiths 
region and the conventional absorbing phase. According 



to eq. (|22|l . the density should asymptotically decay like a 
stretched exponential. To test this behavior, we plot the 
logarithm of the density as a function of t d /( d +^) where 
d = 1 and z « 1.581 is the dynamical exponent of the 
clean one-dimensional contact process. Figure shows 
the resulting graphs for system size L = 10 7 , c = 0.2, 
and several impurity concentrations p = 0.2 ... 0.7. The 
data are averages over 960 runs, each with a different 
disorder realization. The figure shows that the data fol- 
low a stretched exponential behavior In p = —Et - 3875 
over more than four orders of magnitude in p, in good 
agreement with eq. (|22H . The decay constant E, i.e., the 
slope of these curves, increases with increasing impurity 
concentration p. The inset of figure [31 shows the relation 
between E and p = — ln(l —p). In good approximation, 
the values follow the power law E ~ f> z /( d + z ) = ^50.6125 
predicted in 

We now turn to the behavior inside the Griffiths region, 
\° c < A < A c . Figure^shows a double- logarithmic plot of 
the density time evolution for birth rates A = 3.5 . . . 5.1 
and p = 0.3, c = 0.2. The system sizes are between 
10 6 and 10 7 lattice sites, and we have averaged over 480 
disorder realizations. For all birth rates A shown, the 
long-time decay of the density asymptotically follows a 
power-law, as predicted in eq. I|24|) , over several orders of 
magnitude in p (except for the largest A where we could 
observe the power law only over a smaller range in p be- 
cause the decay is too slow). 

The nonuniversal dynamical exponent z' can be ob- 
tained by fitting the long-time asymptotics of the curves 
in figure 01 to eq. I|24l) . The inset of figure 01 shows z' as a 
function of the birth rate A. As discussed in section Hi Dl 
z' increases with increasing A throughout the Griffiths re- 
gion with an apparent divergence around A « 5.2. Unfor- 
tunately, our data did not allow us to make quantitative 



6 



X 


4.3 


V 


5.1 + 


4.1 


o 


4.9 x 


3.9 




4.7 x 


3.7 


A 


4.5 x 


3.5 


□ 




100000 



125 



100 



75 



50 



25 



X 

(top to bottom) 

5.00 

5.10 

5.15 

5.20 

5.22 

5.24 • 

5.26 

5.28 

5.30 

5.32 

5.35 



10 12 
in t 



14 



16 



20 



FIG. 4: Log-log plot of the density time evolution in the Grif- 
fiths region for systems with p — 0.3, c = 0.2 and several birth 
rates A. The system sizes are 10 7 sites for A = 3.5, 3.7 and 
10 6 sites for the other A values. The straight lines are fits 
to the power law p(i) ~ t _1 ^ z predicted in eq. (1241 . Inset: 
Dynamical exponent z' vs. birth rate A. 



comparisons with the predictions for the A-dependence of 
z' because we could not reliably determine z' sufficiently 
close to either the clean critical point or the dirty critical 
point. Close to the clean critical point A°, the crossover 
to the asymptotic power law occurs at very low densities, 
thus the system size limits how close one can get to A°. 
Conversely, for larger A close to the dirty critical point 
A c , the crossover to the asymptotic power law occurs at 
very long times. Thus, the maximum simulation time 
limits how close to the dirty critical point one can still 
extract z' . 



2. Dirty critical point 

After having discussed the Griffiths region, we now 
turn to the most interesting parameter region, the vicin- 
ity of the dirty critical point. In contrast to the clean 
critical birth rate \° c which is well known from the liter- 
ature 43], the dirty critical birth rate A c is not known 
a priori. In order to find A c and at the same time test 
the predictions of the activated scaling picture of section 
III CI we employ the logarithmic time dependence of the 
density, eq. (fT%|) . In figure we plot p^ 1 / 5 with the pre- 
dicted 5 = 0.38197 against ln(t) for a system of 10 4 sites 
with p = 0.3 and c = 0.2. Because the dynamics at the 
dirty critical point is expected to be extremely slow, we 
have simulated up to t = 10 8 (10 9 for the critical curve). 
As before, the data are averages over 480 runs, each with 
a different disorder realization. 

In this type of plot, the logarithmic time dependence 
(|18fl is represented by a straight line. Subcritical data 
should curve upward from the critical straight line, while 
supercritical data should curve downward and eventually 



FIG. 5: p _1//<s vs. ln(t) for a system of 10 4 sites with p = 
0.3 and c = 0.2. The filled circles mark the critical birth 
rate A c = 5.24, and the straight line is a fit of the long-time 
behavior to eq. 11811 



settle to a constant long-time limit. From the data in 
figure we conclude that the dirty critical point indeed 
follows the activated scaling scenario associated with an 
infinite-randomness critical point. The critical birthrate 
is A c = 5.24 ±0.01. At this A, the density follows eq. 1(15)1 
over almost four orders of magnitude in t. 

The statistical error of the plotted average densities 
can be estimated from the standard deviation of p(t) 
between the 480 separate runs. For the critical curve, 
A = 5.24, in figure[Sl the error of the average density re- 
mains below 0.002 which corresponds to about a symbol 
size in the figure (at the long-time end of the plot). We 
have also checked for possible finite-size effects by repeat- 
ing the calculation for a smaller system size of L = 10 3 . 
Within the statistical error the results for L = 10 3 and 
L = 10 4 are identical, from which we conclude that our 
data are not influenced by finite size effects. 



3. Universality 

We now turn to the questions of universality: Is the 
activated scaling scenario valid for all impurity concen- 
trations p and strengths c, and is the value of the criti- 
cal exponent 5 the same for all cases? To answer these 
questions we have repeated the above critical point anal- 
ysis for different sets of the disorder parameters p and 
c. In this subsection, we first show that all these data 
are in agreement with the activated scaling scenario with 
an universal exponent S. We then discuss whether they 
could interpreted in a conventional power-law scaling sce- 
nario as well. 

In the first set of calculations we have kept the impurity 
concentration at p = 0.3, but we have varied their relative 
strength from c = 0.2 to 0.8. In all cases, the density 
decay at the respective dirty critical birth rate A c follows 
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FIG. 6: Time evolution of the density for systems of 10 4 sites 
with p = 0.3 and c = 0.2, 0.4, 0.6 and 0.8 at their respective 
critical points plotted as p -1 ^ vs. ln(t) as in figure^] 

the logarithmic law (|18fl with the predicted S — 0.38197 
over several orders of magnitude in t. Figure shows 
these critical curves for systems with c = 0.2, 0.4, 0.6 and 
0.8. The seemingly larger fluctuations for the curves with 
higher c are caused by the way the data are plotted: The 
large negative exponent — 1/8 strongly stretches the low- 
density part of the ordinate. 

We have obtained analogous results from the second 
set of runs where we kept c = 0.2 constant but varied 
the impurity concentration from p — 0.2 to 0.5. From 
these simulation results we conclude that the data are in 
agreement with the activated scaling scenario for all stud- 
ied parameter values including the case of weak disorder. 
(Note that for p = 0.3, c = 0.8, the disorder-induced shift 
of the critical birthrate is small, (A c — A°)/A° ~ 0.07.) 
Moreover, our results are compatible with a universal 
value of 0.38197 for the exponent 8. 

However, we would like to emphasize that while our 
data do not show any indication of non-universality, we 
cannot exclude some variation of 5 with the disorder 
strength. This is caused by the fact that even though 
we observe the logarithmic time dependence (|18|) over 
almost four orders of magnitude in t, this corresponds 
only to about a factor of 2 to 3 in ln(t). This is a very 
small range for extracting the exponent of the power- 
law relation between p and \n(t). More specifically, the 
asymptotic critical time dependence of the density for 
p = 0.3, c = 0.2 can be fitted by p = {A*\n(t)+B)-°- 38197 
with A m 5.12 and B w —20.1 (this is the straight line 
in figure [SJ). Comparison with eq. I|18|) shows that B rep- 
resents a correction to scaling. For a reliable extraction 
of the exponent one would want the ln(f ) term to be at 
least one order of magnitude larger than B at the very 
minimum. This corresponds to ln(t) 40 or t ps 10 17 
which is clearly unreachable in a Monte Carlo simulation 
for the foreseeable future. 

We now turn to the question of whether the numerical 
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FIG. 7: Replot of the data of figure in log- log form. The 
dashed line shows a power-law decay with the clean critical 
exponent S — 0.159. 



data could also be interpreted in terms of conventional 
power-law scaling. To this end, we first replot the critical 
density decay curves as identified above [44j in standard 
log-log form in Fig. [7] This figure shows that the density 
decay at early times follows a power law with the clean 
decay exponent 8 = 0.159. However, after some disorder 
dependent crossover time t x the curves show pronounced 
upward curvature. This curvature persists to the largest 
times and signifies an asymptotic decay that is slower 
than any power law. The crossover time t x which in- 
creases with decreasing disorder c — > 1 nicely agrees with 
the onset of the asymptotic logarithmic behavior in Fig. 
[B]for all disorder strengths. Thus, the time dependence 
of the density p(t) directly crosses over from the clean 
critical (power law) behavior at short times to the acti- 
vated (logarithmic) behavior (|18f) at large times. We find 
no indication of any power laws other than the clean one, 
not even at transient times. We thus conclude that our 
data, at least for times up to t = 10 s are not compatible 
with power-law scaling at the dirty critical point. 

It would be highly desirable to discriminate between 
power-law and activated scaling by performing a full 
scaling analysis of the data using relations 181911 II) or 
(|15ll(jll7|l . respectively. Unfortunately, such an analysis 
is impossible because the extremely slow dynamics leads 
to the abovementioned strong corrections to scaling in 
the accessible range of times up to t = 10 9 . 



C. Time evolution starting from a single particle 

In addition to the simulations staring from a full lat- 
tice, we have also performed simulations of the time evo- 
lution starting from a single active site in an otherwise 
empty lattice. In these calculations, we have monitored 
the survival probability P s and the number of sites N in 
the active cluster as functions of time. According to eqs. 
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FIG. 8: Survival probability P a and particle number N in the 
active cluster for a simulation starting from single active site. 
The disorder parameters are p = 0.3 and c = 0.2. 



()18|l and (15) , these quantities are expected to behave 
as P s (t) ~ [Init)]- 1 ^ and AT(f) - [ln(i)] e at the dirty 
critical point. 

Figure El shows plots of P~ s and N 1 /® vs. ln(t) for 
a system with p — 0.3 and c = 0.2. The data are av- 
erages over 480 different disorder realizations. For each 
realization we have performed 200 runs starting from a 
single active site at a random position. The system size 
L = 10 6 was several orders of magnitude larger than 
the largest active cluster, thus we have effectively sim- 
ulated an infinite-size system. Figure OH shows that P s 
and N indeed follow the predicted logarithmic laws with 
5 = 0.38197 and 6 = 1.236 at the critical birthrate 
A c = 5.24. Subcritical and supercritical data curve 
away from the critical straight lines as expected. Thus, 
the simulations starting from a single site also confirm 
the activated scaling scenario resulting from an infinite- 
randomness critical point. 



D. Steady state 



FIG. 9: Time evolution of the density in the active phase 
for a system with p = 0.3 and c = 0.2. The system sizes 
are 10 4 for A > 5.8 and 10 3 otherwise. Inset: Steady state 
density vs. distance from critical point. The solid line is a fit 
of the data for A — A c < 0.3 to the expected dirty power law 
Pstot ~ (A — X c f with (3 = 0.38197. The dashed line is a fit 
of the data for A — A c > 0.3 to a power law with the clean 
(3 = 0.2765. 



Therefore, our steady state density values for these A are 
only rough estimates. 

The inset of figure |5] shows the resulting steady state 
densities plotted vs. the distance from the critical point. 
This graph nicely demonstrates the crossover from clean 
to dirty critical behavior. The data away from the critical 
point (A — A c > 0.3) can be well fitted by a power law 
with the clean value of the critical exponent (3 — 0.2765. 
When approaching the critical point the curve becomes 
steeper, and for A — A c < 0.3 the data can be reasonably 
well fitted by a power law with the expected dirty critical 
exponent (3 = 0.38197. The position of the crossover can 
be compared to an estimate from the time dependence at 
A c in figure The critical curve reaches its asymptotic 
logarithmic form at about t x rs 10 4 which corresponds to 
a crossover density of p x sa 0.3 in rough agreement with 
the crossover density in the inset of figure 03 



Lastly, we have studied the behavior of the steady state 
density in the active phase. Close to the critical point 
it is expected to vary as p s t&t ~ (A — Ac)' 5 with (3 = 
0.38197. These calculations require a particularly high 
numerical effort because the approach to the steady is 
logarithmically slow close to the critical point. Therefore, 
we have simulated up to t — 10 9 for birthrates close to 
A c . Figure OH shows the time evolution of the density 
in the active phase for a system with p — 0.3 and c = 
0.2, averaged over 480 disorder realizations (one run per 
realization). The steady state densities can be obtained 
by averaging the constant part of each of the curves. The 
data clearly demonstrate the extremely slow approach to 
the steady state. For A < 5.4 the steady state is not 
reached (within the statistical error) even after t = 10 9 . 



IV. CONCLUSIONS 

To summarize, we have presented the results of large- 
scale Monte-Carlo simulations of a one-dimensional con- 
tact process with quenched spatial disorder for large sys- 
tems of up to L — 10 7 sites and very long times up to 
t = 10 9 . These simulations show that the critical behav- 
ior at the nonequilibrium phase transition is controlled by 
an infinite-randomness fixed point with activated scaling 
and ultraslow dynamics, as predicted in Ref. jsj. More- 
over, the simulations provide evidence that this behavior 
is universal. The logarithmically slow time dependencies 
(|18I19|) are valid (with the same exponent values) for all 
parameter sets investigated including the weak-disorder 
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case. However, the approach to this universal asymptotic 
behavior is extremely slow, with crossover times of the 
order of 10 4 or larger. Since most of the earlier Monte- 
Carlo simulations did not exceed t = 10 5 . . . 10 6 this may 
explain why nonuniversal (effective) exponents were seen 
in previous work. We have also presented results for the 
Griffiths region between the clean and the dirty critical 
points. Here, we have found power-law dynamical behav- 
ior with continuously varying exponents, in agreement 
with theoretical predictions [23, |23 • 

We now discuss the relation of our results to a more 
general theory of rare region effects at phase transitions 
with quenched disorder. In Ref. a general classifica- 
tion of phase transitions in quenched disordered systems 
with short-range interactions 0] has been suggested, 
based on the effective dimensionality d c g of the rare re- 
gions. Three cases can be distinguished. 

(i) If d c g is below the lower critical dimension d~ of the 
problem, the rare region effects are exponentially small 
because the probability of a rare region decreases expo- 
nentially with its volume but the contribution of each re- 
gion to observables increases only as a power law. In this 
case, the critical point is of conventional power-law type. 
Examples in this class include, e.g., the classical equilib- 
rium Ising transition with point defects where d e g = 
and d~ = 1. 

(ii) In the second class, with d e g — d~ , the Griffiths 
effects are of power-law type because the exponentially 
rarity of the rare regions in L r is overcome by an exponen- 
tial increase of each region's contribution. In this class, 
the critical point is controlled by an infinite-randomness 
fixed point with activated scaling. Examples include the 
quantum phase transition in the random transverse field 
Ising model (d e g = d~ = 1) [33, as well as the tran- 



sition discussed in this paper (where d e g — d c = 0). 

(iii) Finally, for d c g > d~ , the rare regions can undergo 
the phase transition independently from the bulk system. 
This leads to a destruction of the sharp phase transition 
by smearing. This behavior occurs at the equilibrium 
Ising transition with plane defects where c? c fr = 2, d~ = 1 
toj, the quantum phase transition in itinerant magnets 
1511 , and for the contact process with extended (line or 
plane) defects [39l loll] . 

Thus, the results of this paper do fit into the general 
rare-region based classification scheme of phase transi- 
tions with quenched disorder and short-range interac- 
tions 45] . These arguments also suggest that the be- 
havior of the phase transition in a higher-dimensional 
disordered contact process should be controlled by an 
infinite-randomness fixed point as well. 

We conclude by pointing out that the unconventional 
behavior found in this paper may explain the striking 
absence of directed percolation scaling [2(J in at least 
some of the experiments. However, the extremely slow 
dynamics will prove to be a challenge for the verification 
of the activated scaling scenario not just in simulations 
but also in experiments. 

Acknowledgements 

This work has been supported in part by the NSF 
under grant nos. DMR-0339147 and PHY99-07949 as 
well as by the University of Missouri Research Board. 
Thomas Vojta is a Cottrell Scholar of Research Corpora- 
tion. We are also grateful for the hospitality of the Aspen 
Center for Physics and the Kavli Institute for Theoreti- 
cal Physics, Santa Barbara during the early stages of this 
work. 



[1] B. Schmittmann and R.K.P. Zia, in Phase transitions 
and critical phenomena, edited by C. Domb and J.L. 
Lebowitz, Vol. 17, (Academic, New York 1995) 

[2] J. Marro and R. Dickman, Nonequilibrium Phase Tran- 
sitions in Lattice Models (Cambridge University Press, 
Cambridge, England, 1996). 

[3] R. Dickman, in Nonequilibrium Statistical Dynamics in 
One Dimension, edited by V. Privman (Cambridge Uni- 
versity Press, Cambridge, England, 1997). 

[4] B. Chopard and M. Droz, Cellular Automaton Modeling 
of Physical Systems (Cambridge University Press, Cam- 
bridge, England, 1998). 

[5] H. Hinrichsen, Adv. Phys. 49, 815 (2000). 

[6] G. Odor, Rev. Mod. Phys. 76, 663 (2004). 

[7] U.C. Tauber, M. Howard, and B.P. Vollmayr-Lee, J. 
Phys. A 38, R79 (2005). 

[8] P. Grassberger and A. de la Torre, Ann. Phys. (NY) 122, 
373 (1979). 

[9] H.K. Janssen, Z. Phys. B 42, 151 (1981); P. Grassberger, 
Z. Phys. B 47, 365 (1982). 
[10] T.E. Harris, Ann. Prob. 2, 969 (1974). 
[11] R.M. Ziff, E. Gulari, and Y. Barshad, Phys. Rev. Lett. 



56, 2553 (1986). 
[12] L.H. Tang and H. Leschhorn, Phys. Rev. A 45 R8309 
(1992). 

[13] Y. Pomeau, Physica D 23, 3 (1986). 
[14] H. Takayasu and A.Y. Tretyakov, Phys. Rev. Lett. 68, 
3060 (1992). 

[15] D. Zhong and D. Ben Avraham, Phys. Lett. A 209, 333 

(1995) . 

[16] J. Cardy and U. Tauber, Phys. Rev. Lett. 77, 4780 

(1996) . 

[17] M.H. Kim and H. Park, Phys. Rev. Lett. 73, 2579 (1994). 

[18] N. Menyhard, J. Phys. A 27, 6139 (1994). 

[19] H. Hinrichsen, Phys. Rev. E 55, 219 (1997). 

[20] H. Hinrichsen, Braz. J. Phys. 30, 69 (2000). 

[21] P. Rupp, R. Richter, and I. Rehberg, Phys. Rev. E 67, 

036209 (2003). 
[22] A. B. Harris, J. Phys. C 7, 1671 (1974). 
[23] A.J. Noest, Phys. Rev. Lett. 57, 90 (1986). 
[24] H.K. Janssen, Phys. Rev. E 55, 6253 (1997). 
[25] A.G. Moreira and R. Dickman, Phys. Rev. E 54, R3090 

(1996); R. Dickman and A.G. Moreira, ibid. 57, 1263 

(1998). 



10 



[26] R.B. Griffiths, Phys. Rev. Lett. 23, 17 (1969). 

[27] A.J. Noest Phys. Rev. B 38, 2715 (1988). 

[28] M. Bramson, R. Durrett, and R.H. Schonmann, Ann. 

Prob. 19, 960 (1991). 
[29] I. Webman et al., Phil. Mag. B 77, 1401 (1998). 
[30] R. Cafiero, A. Gabrielli, and M.A. Munoz, Phys. Rev. E 

57, 5060 (1998). 
[31] J. Hooyberghs, F. Igloi, C. Vanderzande, Phys. Rev. Lett. 

90, 100601 (2003); Phys. Rev. E 69, 066140 (2004). 
[32] F.C. Alcaraz, Ann. Phys. (NY), 230, 250 (1994). 
[33] S.K. Ma, C. Dasgupta, and C.-K. Hu, Phys. Rev. Lett. 

43, 1434 (1979). 
[34] At more general absorbing state transitions, e.g., with 

several absorbing states, the survival probability scales 

with an exponent /3' which may be different from f3 (see, 

e.g., i). 

[35] This relation relies on hyperscaling; it is only valid below 
the upper critical dimension , which is four for directed 
percolation. 

[36] I. Jensen, J. Phys. A 32, 5233 (1999). 

[37] D.S. Fisher, Phys. Rev. Lett. 69, 534 (1992); Phys. Rev. 
B 51, 6411 (1995). 

[38] A.J. Bray, Phys. Rev. Lett. 60, 720 (1998). 

[39] M. Dickison and T. Vojta, J. Phys. A 38, 1999 (2005). 

[40] M.N. Barber, in Phase transitions and critical phenom- 
ena, edited by C. Domb and J.L. Lebowitz (Academic, 
London, 1983), Vol. 8. 

[41] R.S. Schonmann, J. Stat. Phys. 41, 445 (1985). 

[42] R. Dickman, Phys. Rev. E 60, R2441 (1999). 



[43] I. Jensen and R. Dickman, J. Stat. Phys. 71, 89 (1993). 

[44] Identifying the correct A c within the power-law scenario 
is difficult because power laws in pit) are found through- 
out the entire Griffiths region |23l |25| . Our values are 
compatible with the criterion that A c is the smallest A 
supporting a growing N(t) |25|. Moreover, it appears 
highly unlikely that there is a noncritical time scale 
t s > 10 8 after which the curves in Fig. |S| change the 
qualitative behavior they have followed over more than 
four orders of magnitude i n time. 

[45] T. Vojta and J. Schmalian, cond-mat/0405609 to appear 
in Phys. Rev. B . 

[46] In the case of long-range interactions, rare region effects 
can be enhanced further by the interactions between the 
rare regions. This happens, e.g., in the important case of 
magnetic transitions in metals due to the RKKY inter- 
action |47|. 

[47] V. Dobrosavljevic and E. Miranda, Phys. Rev. Lett. 94, 
187203 (2005). 

[48] At a zero temperature quantum phase transition, imag- 
inary time acts as an extra dimension. Because the 
quenched disorder is time-independent, the rare regions 
are effectively one-dimensional. 

[49] R. Sknepnek and T. Vojta, Phys. Rev. B 69, 174410 
(2004); T. Vojta, J. Phys. A 36, 10921 (2003). 

[50] T. Vojta, Phys. Rev. Lett. 90 107202 (2003). 

[51] T. Vojta, Phys. Rev. E 70, 026108 (2004). 



